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We obtain stresses for Newtonian viscous flow in a simple geometry (e.g. 
corner flow, bending flow) in order to study the effect of imposed velocity boun- 
dary conditions. Stress for a delta function velocity boundary condition decays 

as for a step function velocity, stress goes as for a discontinuity in cur- 
vature, tbe stress singularity is logarithmic. For corner flow, which has a 
discontinuity of velocity at a certain point, the corresponding stress has a ~ 

singularity. However, for a more realistic circular-slab model, the stress singu- 
larity becomes logarithmic. Thus the stress distribution is very sensitive to the 
boundary conditions, and in evaluating the applicability of viscous models of 
trench topography it is essential to use realistic geometries. 

Topography and seismicity data from northern Hoshu, Japan, were used to 
construct a finite element model, with flow assumed tangent to the top of the 
grid, for both Newtonian and non-Newtonian flow (power law 3 rheology). Normal 
stresses at the top of the grid are compared to the observed trench topography 
and gravity anomalies. There is poor agreement. Purely viscous models of sub- 
ducting slabs with specified velocity boundary conditions do not predict normal 
stress patterns compatible with observed topography and gravity. Elasticity and 
plasticity appear to be important for the subduction process. 
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Introduction 

The most remarkable and consistent oceanic bathymetric and geoid 
features are found near trenches at subduction zones. One of these features is 
an increase in elevation as the trench is approached from the seaward side. This 
topographic high, or outer rise, is never more than 500 m above the mean level 
of the seafloor. Outer rises start between 1000 km and 3000 km seaward of, and 
crest at about 100 km from, the trench axis. The trench itself is in general a 
narrow low depression about 100 to 200 km in width, 1.5 to 5 km in depth, and 
up to 5000 km in length (Grellet & Dubois 19B2). The landward side of the trench 
is marked by a steep increase in elevation which rises to a crest about 2 to 4.5 
km above the mean level of the seafloor and consists of either continental crust 
or an island arc (Ross & Shor 1965, Hayes & Ewing 1970). Gravity observations 
show variations similar to bathymetry across the trench (Hayes & Ewing 1970, 
Watts & Talwani 1974, Chapman & Talwani 19B2). Other principal observations of 
the features of subduction zones are their geographic distribution over the 
world and the seismicity associated with them (e.g., Uyeda & Kanamori 1979, 
Ruff & Kanamori I960). These features are similar and coherently distributed 
over large scales, and they imply a common cause (Le Pichon et al. 1973). 

Trench-island arc topography, together with seismic Wadati-Benioff zones, 
show that the oceanic lithosphere deforms in a continuous manner as it sub- 
ducts at an oceanic trench. However, understanding the mechanics and dynam- 
ics of the subduction process and their relationship to the topography of sub- 
duction zones and geometry of the Wadati-Benioff zone has proven to be difficult 
(Forsyth, 1980). 

The rheology of the lithosphere and the physical processes associated with 
deep-sea trenches are not yet resolved. For example, is the oceanic lithosphere 
primarily viscous or elastic on the time scale of millions of years ? Various 
theoretical models have been presented to explain the typical shape of plates 
before their subduction at trenches and the distribution of earthquakes 
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associated with subduction in terms of bending. Other than the physical 
features in subduction zones, the surface topography and gravity caused by long 
term (10 e years or more) surface loads such as in seamount chains, isolated 
seamounts, oceanic islands, sediment basins have also been studied using vari- 
ous models. Of these features, the island arc-trench system topography and 
gravity anomalies are significant for understanding the long term and large 
scale behavior of the earth because of their universal distribution and similari- 
ties and their relation to the mantle convection. 

In these models, a mechanically strong lithosphere rests on a much weaker 
asthenosphere. This mechanically strong lithosphere has been modeled as a 
uniform elastic plate (Hanks 1971, Caldwell et al. 1976, Parsons and Molnar 
1976); a uniform plate with elastic-perfectly plastic rheology. (Turcotte et al. 
1978, McAdoo et al. 1978, Chappie and Forsyth 1979); or a uniform viscous plate 
(DeBremaeker 1977, 1980, McKenzie 1977, Melosh 1978. Melosh & Raefsky 1980). 
In addition, various rheological models have been suggested by studies of the 
global gravity field, lithospheric flexure caused by long-term (>10 6 years) surface 
loads, and laboratory experiments (Murrell 1976, Forsyth I960, Watts et al. 1980, 
McNutt 1980). All of these models produce trench topography qualitatively simi- 
lar to that observed. These different rheological models were also proposed to 
explain other field observations and experimental results (Galand 1979, Liu 1980, 
Carey 8c Dubois 1981, Davies 1981). Among these different rheological models, 
viscous models stem from observations of the rate of postglacial rebound, plate 
motion, and evidence of thermal convection in the mantle. (Gutenberg 1958, 
Peltier and Andrews 1976, Isacks et al. 1968, Turcotte and Oxburgh 1967). These 
long term ( >10 8 years) and large spatial scale ( > 1000 km ) phenomena show 
viscous behavior of the earth’s mantle. On the other hand, for short wavelength 
observations, many authors suggest that surface deformation and the gravity 
anomalies near trenches are controlled by elastic forces within the deformed 
plate (Hanks 1971, Le Pichon et al. 1973, Watts and Talwani 1974). 
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However, numerical experiments for mantle convection (McKenzie 1977) 
showed that the surface deformation associated with convection in a fluid with a 
temperature-dependent viscosity bears resemblance to the shape of the ocean 
floor near trenches. DeBremaeker (1977) showed that a viscous model can be 
adjusted to fit the general short-wavelength trench profile fairly well. In his 
model, the shape of a uniform viscous plate is described by a damped sinusoid. 
Results showed that bending stresses are low, and no regional compressive 
stress is required. Melosh and Raefsky computed the stresses in a subducting 
slab using a finite element technique, assuming a viscous constitutive relation 
(Melosh 1978, Melosh & Raefsky 1980). In this model, the lithosphere was divided 
conceptually into two parts, a thin elastic lithosphere and a viscous lithosphere. 
Their principal conclusion is that the first order topographic features of a sub- 
duction zone can be explained by viscous stresses generated during subduction 
in the lower lithosphere beneath a thin elastic upper lithosphere. On the other 
hand, the thin elastic lithosphere has small effect on these features. Other 
models based on a viscous half space were used to study flow-induced pressure 
(Tovish et al. 1978), the angle of subduction (Stevenson and Turner 1977), and 
the compensation of geoid anomalies due to the subduction of oceanic litho- 
sphere (McAdoo 1982). Sleep (1981) developed a numerical model of fluid 
dynamics showing that the gravity field observed over the Aleutians near Adak 
and Amchitka can be explained by a model with power law rheology for the litho- 
sphere near trenches. The underlying assumption in these models is that the 
oceanic lithosphere near trenches behaves as a viscous fluid over the length 
scale of a few hundred kilometers (the scale a < 400 km used by McAdoo 1982) 
and the time scale of a few million years. 

Seismicity in subduction zones (e.g. Ruff and Kanamori 1980) implies that 
the strength of coupling between two converging plates (as measured by M w ) is 
correlated with the age of the subducting oceanic lithosphere and the plate con- 
vergence rate. This may imply that the stress in subduction zones is propor- 
tional to strain rate, which in turn suggests a viscous rheological behavior. 
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In order to better understand viscous flow in representative geometries 
with imposed velocity boundary conditions, we obtain a number of analytical and 
numerical results for simple test problems. We test the appropriateness of the 
viscous lithosphere model by investigating some simple test problems as well as 
more realistic models. The reasons to do this include (1) no study has used 
actual topographic observations to develop theoretical models and (2) no one 
has demonstrated that the results obtained from any particular theoretical 
model are comparable with the realistic case. Later we will show that the stress 
distribution in viscous flow is very sensitive to velocity and stress boundary con- 
ditions. This indicates that it is very important to use observations near subduc- 
tion zones to construct a viscous flow model for a given individual subducting 
slab rather than using general theoretical expressions. In order to study the 
stress distribution in a flow model of realistic geometry of a subducting slab, we 
use the finite element method, since it has advantages over other methods for 
models with arbitrary geometry. 

We use the topography and seismicity data from the Japan trench to pro- 
vide constraints on the geometry of finite element models. The Japan trench is 
one of the best studied trenches. The bathymetric data obtained during 1964 - 
1980 (Tomoda and Fujimoto 1982) were used to get an average profile transverse 
to the Japan Trench near northern Honshu. The linear Wadati-Benioff zone in this 
region was used to determine the location of the subducting lithosphere. The 
gravity data in this region clearly resemble the topography data. A detailed 
analysis of focal mechanism solutions of small and intermediate size earth- 
quakes, together with the seismic coupling (measured by M w ), also contributed 
to our understanding of the stress state in the slab and the rheological proper- 
ties of the upper mantle in subduction zones (Kawakatsu & Seno 1983, Ruff and 
Kanamori 1980, Kanamori 1971, 1977). 

We find that the purely viscous models of a subducting slab do not predict a 
stress pattern compatible with the observed topography and gravity anomalies 
over the Japan trench. The results from simple test problems indicate that the 
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discrepancies between observations and numerical results for modeling sub- 
ducting slabs are general for viscous models. Based on these calculations, we 
conclude that purely viscous models for the lithosphere do not predict a stress 
pattern consistent with the actual observations of topography and gravity. 

In the following sections, we calculate pressure distributions for the flow 
induced in a viscous incompressible fluid by the motion of boundaries. First, we 
discuss half space flow because its simple geometry, then calculate numerical 
models for more realistic geometries. 


Analytic formulations 

To understand viscous flow with imposed velocity boundary conditions, we 
first discuss some analytic results. Later on, we use these results to qualita- 
tively explain the numerical results. 'Hie equations of motion for a uniform den- 
sity, incompressible fluid undergoing steady flow are 


( 1 ) 

V • ? = 0 (2) 

where r) is the viscosity, p the density. ^ the fluid velocity, and p the nonhydrosr 
tatic or hydrodynamic pressure. In a homogeneous medium, both the density 
and the viscosity are taken to be constant. If the flow is assumed to be two- 
dimensional in the x.y plane, it can be described by a stream function, V*. 




dji _ 9 J 
9y' 8x’ j 


where = V'(x.y) satisfies the biharmonic equation 


(3) 


vv = o. 


(4) 


Any solution of (4) must have the form 


^ = d + 0x + y y 


(5) 



■where a . /3 and 7 are harmonic (Pearson 1959). For our problems, it is con- 
venient to use the following representation for the stream function: 


^(x,y) 


= |a si 


sinh ky + B cosh ky + Cy sinh ky + Dy cosh ky 


sin kx 
cos kx 


( 6 ) 


where A.B.C.D are constants to be determined, and k is a wavenumber. Using (3) 
and (l) both the velocity and the pressure associated with the flow can be deter- 
mined. Other stress components may then be obtained by : 


: = ~P + 2r) 


fix 


(?a) 




(7b) 


°yy = -p + 27? 



(7c) 


The total stress may be obtained by summing the hydrodynamic stress and the 
hydrostatic stress. The velocity and hydrodynamic stress components 
corresponding to the flow described by the stream function (6) are given as fol- 
lows : 


v,(x.y) = |E sinh ky + F cosh ky + Dky sinh ky + Cky cosh ky 


sin kx 
cos kx 


( 8 ) 
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p(x,y) = p(x.-«) + 2r)k 


C sinh ky + D cosh ky 


— coskx 
I sin kx 


( 10 ) 




= - 2rjk 


sinh ky + H cosh ky + Dky sinh ky + Cky cosh ky 


-cos kx 
sin kx , 


-p(x,-») 


(ID 


a *s - 2t7k|r sinh ky + E cosh ky +Cky sinh ky + Dky cosh ky 

°yy = ~2p ~Osx 

where E = C + Bk. F = D + Ak. G = 2C + Bk, H = 2D + Ak. 

Two dimensional flow in a half space 

We consider two-dimensional flow in the half space y £ 0, — <»<x<°°, for vari- 
ous simple velocity boundary conditions applied at the surface, y = 0. We 
assume velocity goes to zero as y approaches — *>. The various velocity boundary 
conditions are: i) a harmonic function; ii) a delta function and; iii) a Heaviside 
unit step function. Using formulae (B) to (13), we obtain some useful results 
summarized in Table 1, which include stream functions and pressures. Other 
stress components may also be obtained by these formulae. The fundamental 
results shown in Table 1 indicate that a point source of horizontal velocity on the 
surface of a half space generates a region of compression in the forward direc- 
tion and a region of tension in the backward direction. On the other hand, a 
point source of upward velocity generates a region of tension beneath this point 
source and a region of compression elsewhere. A horizontal velocity 


(13) 


sin kx 
icos kx 


( 12 ) 
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discontinuity (step function) produces a region of tension below the discon- 
tinuity since it creates a divergence in surface velocity. A discontinuity in the 
vertical velocity produces an antisymmetric pressure distribution, which is ten- 
sional beneath the region with outward velocity. 

Although the geometry of the flow is very simple, we found that the results 
obtained are conceptually useful for understanding the stress distribution of 
flow with an arbitrary geometry calculated using finite element methods. In 
order to see how the stress changes due to the variation of the geometry of flow, 
we study corner flow next for comparison. 

Corner flow 

The basic corner flow model is discussed by many authors (Batchelor 1970, 
McKenzie 1969, Stevenson and Turner 1977, Tovish et al. 197B). We assume that 
the velocities on the surface of the flow’ are given, and use the same method to 
calculate the pressure within the corner. The geometry of the flow is shown in 
Figure la. The velocity vectors on the surface have both tangential and normal 
components. In polar coordinates the stream function is given by 

’4>{r,0) - r (A sin 6 + B cos 8 + CS sin 8+ D£cos 0), (14) 

where the constants A, B, C, D are to be determined by boundary conditions. 
Assume that the boundary conditions for the flow within region A (Figure la) are 


^ = A 3 , + a 3 $ 

on 8 - 0 

(15a) 

$ = ft 3 ,. + (} 

on 8= 8m 

(15b) 


where ^ is a radial unit vector and 3$ the azimuthal unit vector. The velocity 
may be written as 

(16) 


The stream function (14) satisfying the boundary conditions has constants A - D 
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given by 


A = 


a6 m - f}8 m cos8 m + asin^cos^ - ^sinfl m - X 0 Z - (j.8 m sin6 a 


(17a) 



(17b) 


C = 


£8 m sin8 m - as'm 2 8 m - X(sin0 m cos0 m -0 m ) + ^(sin# m - 8 m cos6 m ) 


(17c) 


D = 


1 _ 

A 


- a8 m + 06, 


m cos 8 m - asin£? m cos0 m + /?sin 6 m + Xsin z 0 m + fi8 t 


- m sin 0 m r 


(17d) 


A = sin 8 0 m - 8^ 


(I7e) 


For a Newtonian fluid, the stream function can be written as t(/= r0, and the 
pressure is given by 

p = — n_( 0 '“ + 0 ) = ?2L(c sin 8+ D cos 6) (18) 

r r 

(Tovish et al. 1978). Thus, the pressure corresponding to the boundary condi- 
tions (15) is 

P = EPi (1 9 ) 

t=i 


2tt_ X 

Pl= rT 


6 k,sin 0 + sin£ m sin( 0 m - 6 ) 


(19a) 


w r A 


sinf? m sin£? + $nSin(d? m -e) 


(19b) 
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(19c) 
(I9d) 

Since (15) represents any velocity discontinuity for corner flow, the pres- 
sure (19) singularity for a velocity discontinuity is of order l/r. The results 
agree with those of Table 1 for 6^ = tt. 

A circular-arc slab model 

A more realistic model for subducting slabs is the circular-arc slab model. 
The surface of the slab is composed of two straight lines on either side of a cir- 
cular arc of angle 8. This model is shown in Figure lb, where the subducting 
slab is a strip of constant thickness in the x, y plane. Using the finite element 
method, stresses may be obtained for any geometry or boundary conditions. 
However, for comparison with previous results for half space flow and corner 
flow, we present an approximate solution for the circular-arc slab model. 

Assume the viscosity is constant in the slab, and the velocity along the top 
surface, which is composed of two straight lines and an arc, is constant and 
tangential to the surface along its length. Consider the velocity boundary condi- 
tions on the top surface of the slab 

v x = — v; v y = 0 at x &0, y = 0 ( horizontal section ) (20a) 

v, = — vcos <p ; v y = -vsin <p at Vx 2 + (y+b) z =b, xfiO, y<8 ( arc section )(20b) 

v x = —v cos 8 ; v y = — v sin 8 at xfiO, y = 0, p ^0 ( dipping section ) (20c) 

where <p is the angle from the beginning of the bend of the slab ( x = 0 ) and v is 
the net convergence rate in the subduction zone. The top surface of the slab in 
this model has only two changes in curvature. This is a simple model by which 
we can obtain some insight into the relationship between the stress distribution 


Pa 


__ —2.7) a 


r A 


l sin8 m cos(8 m - 8) + 0 m cos( 


P4 


= 2tl i_ 

r A 


sin0 m cos£?+ 8 m cos(8 m -8 ) 
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and the curvature change. Assume the bottom surface of the slab is traction 
free. 

An exact analytical solution of this problem has not been obtained. A 
numerical solution for this problem can be obtained by using finite element 
technique (Melosh & Raefsky 1980). However, by looking at some simple cases, 
for which we have analytic solutions, it is possible to estimate the stress in a 
more complex geometry. For the present problem, it is important to estimate 
stresses near the bending points C: x = 0, y = 0, and D: x = - b sin# , y = b cos#- 
b. The following two stream functions describe the flow near these two points. 
First, the flow in the circular section arc has the stream function 

Y' = Cjr 2 (21) 

and, secondly, the flow in the straight sections has the stream function 

Y' = c 2 r sin(<p + a) (22) 

where Cj,c 2 , and a are constants. The two stream functions do not match at 
points C, D, E, and B, and therefor stress discontinuities appear near these 
points. 

We first show analytically that the singularity in pressure is logarithmic at 
the bending points C and D. To satisfy the three different boundary conditions 
(20), the stream functions (21) and (22) become 

Y'(x.y) = -v(y+b) (23) 

Y'(x.y) = -^-[x 2 + (y+b) 2 ] (24) 


Y(x.y) = —v [(y+b) cos# - x sin#] 


(25) 
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respectively. Substituting these three stream functions into (3), we get the velo- 
city at y =0 

v x = -v; v y = 0 at x ^ 0, y=0 (26a) 

v x = — v, v y = v -- at — b tan x ^ 0, y=0 (26b) 


v x = -v cos# Vy = —v sin# at x< — b tan#, y=0 (26c) 


If the dip angle of the slab is small or we are only concerned with the region 
near bending point C, we may assume (26) as the velocity boundary condition for 
a half space. Using the velocity boundary conditions for the half space, we can 
make use of the results in Table 1, which we take as Green functions. This is 
feasible because the problem is linear. Thus, the pressure p(x,y) ( — y^O 
) can be obtained by superposition: 


( \ 2*7 v r C [y 2 - (x-xp) 2 ] x 0 dx 0 

P1X,y; ' rrb '-bun* jy + (x-xo) 2 ] 


— v (1-cos 6 ) (2-) 

TT 


2y 

(x+btan # ) 2 +y* 


— v (tan 0-sin 9 ) ( 

7 r 


— 2(x+btan #) 
(x+btan 6 ) 2 + y 2 


(27) 


The three terms in (27) are due to the changes of vertical velocity on the seg- 
ment — btan#^x:s0, y=0 defined by (26b), the discontinuities of horizontal 
and vertical velocity at x = -b tan# y=0 respectively. The value of the integral 
is: 


. _ 1_ y 2 + (x+btan#) 2 btan# (x+btan#) 

2 n s^+y 2 (x+btan#) 2 +y 2 

The pressure is 

2 77 V (x+btan#) (2tan#-sin#) — y (1 -cos#) 
P= * (x+btan# 2 + y 2 


( 28 ) 
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ZLX_ ln f. + ( x + , b tan 

tt b x z + y 2 


(29) 


There are other ways to estimate stresses in the slab but the above approxi- 
mation is simple and adequate for estimating the stress in the region near the 
singular point x = 0, y = 0. The approximation does not predict the velocities 

6 6 

defined by (20) on the section of as accurately as those where 

w Ct 

Hence this expression is only an approximate solution for the pressure in the 
region near the singular point x = 0, y = 0, and is ineffective for the region 

Q 

tp The equation (29) shows that the pressure is logarithmically singular at 

the point x = 0, y = 0. Since the problem for circular arc is antisymmetric 

Q 

about tp - the pressure associated with the condition (29) has two loga- 
£ 

rithmic singularities at point C and D on the top surface of the slab. Formula 
(29) shows that the logarithmic singularity is tensional at point C, and compres- 
sional at point D by the antisymmetry of the problem. 

Using finite element methods, pressure profiles along lines near the top sur- 
face for this slab model have been calculated for comparison with the approxi- 
mate solution (29). Figure 2 shows the pressure profiles along two lines parallel 
to the top surface of the slab calculated using formula (29), P a , and the finite 
element method, P f . In figure 2, the distance d between each line and the top 
surface is given in units of the thickness H v of the slab, as is the horizontal dis- 
tance measured from the beginning of the bend, point C. Pressure has units of 
i}\/ Hy. Near point C, the pressure profiles from the finite element calculation P f 
and approximate solution P a agree well. Although the value of the pressure at 
point C for the approximate analytic solution, P a , is infinite, the values outside a 
small neighborhood of point C agree well with Pj. Since we have ignored the 
influence of the lower boundary of the slab, (29) is only valid near the surface of 
the slab. This is why we did not get a singular point at point B. 
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From the stresses obtained in Table 1, those for corner flow and those for 
bending flow, we see that the stress distribution is very sensitive to imposed 
velocity boundary conditions. Stress for an impulsive velocity boundary condi- 


tion behaves as for a step function velocity, stress goes as and for a 

discontinuity in curvature, the stress singularity is logarithmic. For corner flow, 
which has a discontinuity of the vector velocity at certain point, the correspond- 


ing stress has a — singularity at that point. However, for the more realistic 
w r 

circular-slab model, the stress singularity becomes logarithmic. The geometry 
of seismic Wadati-Beniofl zones around the world shows that the curvature of 
subducting slabs probably varies from the outer rise down to the upper mantle 
(Isacks & Barazangi 1977). Since the stress pattern is sensitive to the change in 
curvature, the comer flow model is probably not a good model to study detailed 
stress distributions in subduction zones. 


Finite element calculations for Helosh & Raefsky '8 models 

The problem of subduction of a viscous slab may be solved numerically, as 
was done previously by Melosh Sc Raefsky, 1980. Melosh (1978, Melosh Sc Raefsky, 
1980) assumed that the lithosphere could be treated as two parts : an elastic 
upper part, 15 to 25 km thick (H e ), and a viscous lower part 50 - 80 km thick 
(H v ). The asthenosphere is a zone of low viscosity lying below this lithosphere. 
Melosh and Raefsky assumed that the elastic bending stresses are unimportant, 
and elastic stresses are not included in their calculations. 

Their calculations show the following results. As the lithosphere subducts 
and bends downward, extensional flow in the upper part of the very viscous layer 
produces a zone of extensional normal stress. If, as they suppose, the stress 
propagates through the elastic layer to the surface without significant changes, 
this extensional stress will produce a topographic low, where the oceanic trench 
is located. The topographic high on the landward side of the trench is due to a 
zone of compressive stresses. As their model produced qualitatively reasonable 
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results, we wished to investigate it further using more realistic geometries. We 
first attempted to reproduce their results. 

Melosh and Raefsky (1980) used two finite element grids. One has continu- 
ous curvature along its top surface with error function derivatives, while another 
is the circular-are grid, with the geometry of the circular-arc slab model we 
have used above. Since the numerical results for these models are important 
for modeling the subduction zone, we discuss these models in detail. We studied 
these models by using finer grids than that shown in Figure 2 of their paper, and 
analyzed the results for different boundary conditions. Following Melosh and 
Raefsky, we do not include the elastic layer in our present calculations. While 
they used a Maxwell viscoelastic rheology and calculated its steady-state 
response, we posed the problem as incompressible creeping flow and solved the 
problem using a penalty function formulation (Hughes Sc Taylor 1978). 

The boundary conditions for the viscous part of the lithosphere are as fol- 
lows. The base of the lithosphere is free of stress other than hydrostatic pres- 
sure, and the velocities on the top of the viscous lithosphere are specified. 

To have adequate resolution of stresses for the circular-arc slab model, we 
used a fine grid. The numbers on the Figure 3a give the numbers of elements 
along each part of the grid. A total of 15 x ( 16+8+45+19+10 ) = 1470 bilinear 
isoparametric elements were used. The velocity boundary conditions (20) are 
applied to the grid such that the velocity vectors are tangent to the upper line 
of the grid. The left end of the slab is subject to constant velocity. The right 
side of the slab has the same velocity boundary condition as that of the horizon- 
tal section. These are given as: 

v x - — v, v y = 0 at the right side of the grid (30a) 

v x = —v cos#, Vy = -v sin# at the left side of the grid (30b) 

The results for the calculated pressures are shown in Figure 3. Figure 3a 
shows the pressure contours for the model. The contours are symmetric and 
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have four singularities at points B, C, D, and E (Figure lb). The two singularities 
on the top are logarithmic. Figure 3b is a plot of the pressure along the top line 
of nodes of the grid. 

The derivative of the prescribed velocity (20) is discontinuous at the begin- 
ning and end of the bend. To alleviate this discontinuity, Melosh & Raefsky 
(1980) suggested employing a smoothed velocity profile in the circular arc sec- 
tion given by 

'^|S£.[l-co.(Sg4] 

v a = -Vv 2 — v| (31b) 


0 <<p< 


e 


T 


2 
<p<B 


(31a) 


Using the velocity boundary condition above, the stress curves calculated 
along a line parallel to the top surface of the slab have a small uplifted region or 
bench. Since the bench was not fully resolved in Melosh & Raefsky’s paper 
(1980), we calculated the stress using a finer grid. The grid used has 1400 ele- 
ments. and the size of the elements near the edges of the grid is much smaller 
than that of the elements near the center line. Figure 4a shows the pressure 
contours for the circular-arc grid with smoothed velocity boundary conditions 
(31). Note the differences between the pressure contours in Figure 3a and 4a. 
In Figure ,4a. the maximum and minimum values of pressure move toward the 
bend axis 1 and spread out, loosing their sharp points. Figure 4b is the pressure 
profile along the top layer of the finer grid. This pressure profile is smooth and 
does not have sharp points at the two ends of the bend. By inspecting similar 
plots for nodes at depth in the slab, we find that the amplitude of the bench near 
the bend axis is most pronounced at the surface and decreases with depth. 

In the stress analysis for the corner flow and the circular-arc slab model, as 
shown in previous sections, discontinuities of velocity or its derivatives produce 
local stress variations. The circular-arc slab model has discontinuous 
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derivatives of vertical velocity, and thus its curvature. The stresses for the 
model have logarithmic singularities. In this case, we have different velocity 
profiles applied to the surface of the grid on each side of the bend axis, where 
they have discontinuities of the second derivatives of velocity. These discon- 
tinuities also produce local variations of stress, which causes the bench in the 
stress profiles. 

Comparing the pressure distributions for the circular-arc slab grid with the 
two slight different velocity boundary conditions we used above, we see that for a 
particular grid, the stress is very sensitive to velocity boundary conditions. In 
order to study the effect of different grid geometries, we also investigated a 
finite element model with an error function velocity profile, given as 

v , = -V^_ cr .c|iL] 

(32) 

V K = -V^-Vy 

I 

where S is the arc distance from the bend and the cut-off parameter S c is given 
by 

S c =^-sin(|j (33) 

Using the same boundary conditions on two different grids, we find that 
stress distribution patterns are quite similar. One grid is the circular-arc grid 
used previously (Figure 3a), and the other is a rectangular grid, i.e. a flat slab 
with no dipping region. Figure 5(a-l), 5(a-2), and 5(a-3) show 3 different velocity 
profiles which correspond to 3 different ’cut off’ parameters of the error func- 
tion velocity (32) and (33). Applying these velocity profiles to the circular-arc 
grid in Figure 3. we have pressure profiles at the surface, Figure 5(b-l),5(b-2), 
and 5(b-3). For a rectangular grid, we have pressure profiles 5(c-l), 5(c-2), and 
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5(c-3). Comparing these pressure curves and velocity profiles, we notice that for 
the same velocity boundary conditions, the difference in the magnitude of the 
maximum tension is about 20 percent for the 2 different grids. However, the 
shape of the pressure curves is quite similar. Since, for the same grid, the 
stress patterns are very different for different velocity conditions, the stress 
variation is more strongly controlled by velocity boundary conditions than by 
the geometry. 

According to equation (9), the velocity gradient V $ causes stresses in 
viscous flow. From the results in Table 1, we can understand in a general way 
the results obtained by the finite element method. We consider the velocity vec- 
tor along a particular stream line. At each point along the stream line we have 
two characteristic directions, tangent and normal to the stream line, respec- 
tively. We project the velocity vectors of a small segment of the stream line 
onto these two directions, and call the variations of the two components along 
the streamline tangential and normal gradients. The gradients defined above 
are discrete; for example, the tangential and normal gradients are zero on the 
right side of the slab of point C in Figure lb. The normal gradients are negative, 
the tangential gradient is zero on the entire circular arc for constant tangential 
velocity boundary conditions. At point C and D, we have nonzero normal and 
tangential gradients. From the fundamental solutions for the half space flow in 
Table 1, we see that tangential gradients cause tension, while normal gradients 
cause compression followed by tension. Since the fundamental results referred 
to in Table 1 are for a half space, this kind of intuition is effective only for flow 
like our thick plate model. 

The stress distribution obtained from the finite element calculations, in 
particular, the stress variation near the bending point, can be explained using 
the analytical results. In Figure 3, the pressure is compressional to the right of 
the bending point C. Near C, the pressure changes sign becoming tensional and 
then increases rapidly to a tension maximum at C. 
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Before the bend, the normal and tangential gradients are both zero. After 
the bend begins, the normal gradient of velocity becomes nonzero, being nega- 
tive to the left of point C. But the tangential gradient is negative near C and 
zero to the left of it (Figure 6). This is the reason that the maximum tension 
appears at the bending point C. The compression before the point C is caused by 
the negative tangential gradients near point C. The pressure distribution near 
the end point, D. of the arc can also be explained by negative normal and posi- 
tive tangential gradients. 

Figure 7 shows the superposition of pressure due to the normal velocity 
gradients for a series of vertical step function velocity discontinuities. The rea- 
son that the pressure near the bend axis tends to be zero in the unsmoothed 
boundary condition calculations is that, although the normal gradient there is 
negative (Figure 7), the total pressure will vanish near the bend axis by superpo- 
sition of the stresses due to a series of vertical step function velocity discon- 
tinuities. 

The above discussion can be extended to explain the 'bench' phenomenon 
found in Melosh & Raefsky 's (1980) finite element calculations. From the 
analysis of the circular-arc model, we know that by applying a tangential velo- 
city boundary condition, the model produces reasonable results, which match 
the approximate analytic solution. In Figure 3, the bench does not occur. How- 
ever, Figure 4 shows a bench near the bend axis. The only difference between 
these two models is that slightly different velocity boundary conditions are used. 
Subtracting the two velocity profiles, we obtain a differential velocity boundary 
condition. Applying this boundary condition to the same grid, we obtain the 
results shown in Figure 8(b-2). 

Figure 8a shows three velocity boundary conditions for the circular-arc 
model of Figure 3a. Since the smoothed velocity boundary condition differs 
from the unsmoothed velocity boundary condition only on the circular arc parts 
of the grid, only the velocities of smoothed, unsmoothed, and their difference 
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are plotted. They are unsmoothed vertical velocities A y and horizontal velocities 
A,, smoothed vertical velocities B y and horizontal velocities B x . and their 
differences C y , C x . The differential velocities C are also plotted on coordinates 
with the origin at the bend axis and a 45° rotation of the original coordinates 
(Figure 8(a - 2)). Note that 8(a - 1) shows that the unsmoothed vertical veloci- 
ties Ay and differential velocities C y have derivative discontinuities at the begin- 
ning and the end of the bend which have the same magnitudes but reverse signs. 
Figure 8( b - 1 ). 8( b - 2 ), and 8( b - 3 ) show the pressure profiles along the top 
line of the grid. They correspond to unsmoothed, differential and smoothed 
velocity boundary conditions respectively. Both unsmoothed and differential 
velocity boundary conditions produce logarithmic singularities at the beginning 
and end points of the bend. They have opposite signs and are canceled in the 
calculations with smoothed velocity boundary conditions. The bench occurs in 
the calculations with the differential velocity boundary condition, since it arises 
in calculations with the smoothed boundary condition. 

In order to explain the occurrence of these three pressure profiles intui- 
tively, Figure 8 shows a rough scheme explaining how the smoothed pressure is 
obtained. B( c - 1 ) shows the normal gradient gi of velocity on the arc, which is 
a (negative) constant for the unsmoothed velocity. B( c - 2 ) is the pressure 
profile along the top line of the grid produced by the velocity. 8( d - 1 ) shows 
the residual or differential normal velocity gradients g' a corresponding to the 
differential velocity. 8( d - 2 ) shows the two logarithmic singularities and a 
bench produced by the residual normal gradients. 8( a - 1 ) and B( a - 2 ) also 
show that the tangential velocity gradient is very small compared to the normal 
gradient for the differential velocities, and does not make as significant a contri- 
bution to the stress on the arc section as the normal gradients do (Figure 8c 
and 8d). The important finding is that the ’bench’ is not the result of sharp 
bending and unbending of the slab, but rather is the result of a small kink in 
velocities, resulting from smoothing the velocity boundary condition near the 
axis of the bending slab. The choice of .the smoothing function produces a 
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’bench’ which severely distorts the stress distribution near the bending axis at 
shallow depth. One must therefore be extremely careful in interpreting the 
results of numerical models in a subduction zone type geometry; small 
differences in velocity boundary conditions affect the stress qualitatively as well 
as quantitatively. 

ilodels using observed grids for the Japan trench 

Prom the above discussion, we can see that for the flow models addressed, 
the stress distribution is very sensitive to velocity boundary conditions. Thus, 
instead of using a simple theoretical model, it is necessary to construct viscous 
flow models with more realistic geometries for individual subducting slabs. 

The path of the subducted plate in the upper mantle is usually delineated 
by a zone of earthquakes which occur in the upper part of the plate, but the 
detailed relationship between the upper surface of the descending plate and the 
seismic zone may not always be clear. However, for the Wadati-Benioff zone near 
the Northern Honshu Arc, Japan, detailed seismicity studies show that earth- 
quakes with the same mechanism are systematically distributed (Kawakatsu et 
nl. 1983). 

We assume that the planar structure of seismicity represents the geometry 
of the subducting plate. This is probably true for thrust type earthquakes which 
have a seismic slip front identical to the aseismic front defined as the seaward 
edge of the aseismic mantle wedge beneath the island arc (Yoshii 1975). Thus 
the approximate location and dip of the interface between the two converging 
plates are indicated by the shallow thrust zone and deep thrust zone (Kawakatsu 
and Seno 1983). We assume that there is no sharp bending of the oceanic plate 
in the asthenosphere, and consider the upper envelope of the epicenters of man- 
tle earthquakes as the upper surface of the subducting plate. These constraints 
are similar to those of Isacks and Barazangi (1976). The thrust and compres- 
sional type earthquakes in cross sections of the whole Northern Honshu Arc were 
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digitized. A polynomial curve was generated by a least squares fit to these epi- 
centers to obtain the geometry of the interface between the two converging 
plates. Another polynomial, which fits the averaged topography data, has a 
small oSset 50 km landward from the trench from the polynomial which fits the 
Wadati-Beniofl zone (1.154 km below the topography polynomial). In obtaining 
the least squares polynomial, we did not consider the data points near the 
seamount (Tomoda and Fujimoto 1982). The profiles of topography correspond- 
ing to the boundaries between the four seismic sections of Kawakatsu and Seno 
(1983) were obtained by digitizing the bottom topography map in this region 
(Tomoda and Fujimoto 1982). We construct a finite element grid assuming that 
the thickness of the viscous lithosphere (H v ) is constant and equal to 50 km. 

The boundary conditions for the models are as follows: The base of the 
lithosphere is free of shear stress and normal stress other than that due to 
hydrostatic pressure. At the top surface of the viscous lithosphere, the velocity 
is constant and tangent to the surface of the lithospheric slab. For the down 
going side of the viscous lithosphere, at the depth of the deepest earthquake in 
the data (about 180 km deep), the velocity is parallel to and equal to the velocity 
on the surface of the lithosphere. On the oceanward edge of the viscous litho- 
sphere. the velocity is also defined as a constant and parallel to the upper sur- 
face. 

We assume that the steady state bending of the elastic plate forces the 
upper surface of the viscous layer to move with a constant tangential velocity, 
also assuming that the thickness of the upper elastic layer does not change. We 
note that the double seismic zones are almost parallel over 200 km from the 
trench. This may imply that there is no significant change of the thickness of 
the subducting slab. Although the velocity along the upper surface may change 
from place to place because of coupling and uncoupling of the two converging 
plates, as a first approximation, we assume it is constant and tangent to the sur- 
face, as the basic question here is, to first order, whether or not the topography 
of the island trench-outer rise system can be supported by the viscous stresses 
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in the flow of a viscous layer beneath an elastic lithosphere. The free stress con- 
dition at the base of the viscous layer is appropriate if one considers the 
asthenosphere as a zone of low viscosity underlying the lithosphere. 

Since the position of the interface between the two converging plates is 
uncertain, several grids which have slightly different geometries of their upper 
boundaries were used. These boundaries are determined by polynomials which 
fit the upper envelope of the hypocenters somewhat differently. These grids 
differ slightly from each other with respect to bending curvature and depth of 
the interface of the two converging plates. Finite element calculations show 
that although these models yield somewhat different results, none adequately 
predict the trench topography. As the implications from different models are 
quite consistent, we will discuss those features of the models which do not 
depend on the geometry of the grids. Figure 9 shows two of the models used in 
our finite element calculations with slightly different geometry. The models 
shown in Figure 9a and 9b, will be called model A and B respectively. 

The upper boundary of the grid of model A is an estimate of the top surface 
of the descending lithosphere based upon the following constraints. The bottom 
topography of the sea is taken to be the upper surface of the plate just before 
thrusting beneath the landward plate. Figure 10 shows a bathymetric profile 
east of Northern Honshu as a function of distance from trench. It ranges from 
the trench ( x = 0.0 ) to 900 km ( x = 1B.0 ) seaward. It was obtained by averag- 
ing 5 cross sections of topographic profiles from 450 km west to 900 km east of 
the trench axis. These cross sections correspond to the boundaries of 4 seismi- 
city cross sections used by Kawakatsu and Seno (1983). The top lines of the two 
finite element grids, for model A and model B are also plotted in Figure 10. The 
two curves (upper surface of the grid) fit the topography very well in general. At 
500 km from the trench there is a 50 km wide 300 m high upper swing. It results 
from averaging the topography data including a seamount in the fifth section. 
Except for this, the sea floor is quite flat at that distance. Figure 10 also shows 
the outer rise is 500 m above the sea floor and 500 km wide. Model B differs 
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from model A slightly in curvature. 

Observed grid with constant Newtonian viscosity 

We present the stresses calculated for model A, assuming a constant 
Newtonian viscosity. Figure 11 shows the stress components T^.Tyy, and pres- 
sure P as functions of distance from the trench for the top surface of the model. 
To calculate the predicted topography from the viscous stresses, we assume a 
simplified situation in which the normal stresses Tyy at the top surface of the 
viscous lithosphere are fully compensated by a decrease in elevation of the over- 
lying rock: 

W = -T w /(p m -p w )g (34) 

where W is the decrease in the depth of the sea floor, Tyy is the normal stress at 
the base of the elastic lithosphere. p m is the mantle density, p w is the density of 
the seawater and g is the acceleration of gravity. In this case, the normal stress 
T,y produced by flow in viscous lithosphere is proportional to the depth of the 
sea floor. For the vertical stresses Tyy from BOO km seaward to 250 km land- 
ward of the trench axis obtained in our model, there is a maximum tension of 
7xl0 fl Pa, 200 km landward of the trench and a maximum compression of 
5x10° Pa, 140 km seaward of the trench. Denoting the ^ and x,. as the distances 
from predicted trench axis and outrise to the actual trench axis, W t and W r as 
the corresponding topography from sea floor (positive up) and assuming p m = 
3300 kg m~ 3 , p w = 1000 kg m -3 , g=10ms -z , we obtained the predicted topography 
x t = — 200 km, W v = - 300 m, x T = 140 km, and W r = 20 m. This does not agree 
with the actual topography even qualitatively. 

The above assumption neglects the ability of the elastic part of the litho- 
sphere (the upper 15 - 30 km) to transmit normal forces to moderate distance. 
For an elastic lithosphere loaded by vertical stresses Tyy applied at its base, the 
deflection W can be calculated by: 
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W(x) = ~§D~f-~ e “ ( cos ^“ + sin--— ) Tyy(s) ds 


a 


(35) 
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E is the Young’s modulus, H e is the thickness of the elastic lithosphere, v is 
Poison's ratio, D is the flexural rigidity and a is the flexural parameter (Turcotte 
and Schubert 1982). The elastic part of the lithosphere can be considered as a 
filter with flexural parameter usually assumed to be 40 to 75 km (Walcott 1970, 
Watts et al. 1975, McNutt & Menard 1978). Deflections due to normal forces with 
wavelengths much longer than the flexural parameter are transmitted through 
the elastic lithosphere essentially unaltered. 

Using the viscous stresses obtained in our model, we calculated the 
deflection of the elastic lithosphere by formula (35). The viscous stress Tyy 
obtained from our model has wavelength about 500 km (Figure 11), so it can be 
transmitted through the flexible elastic lithosphere resulting in the surface 
topography proportional to the stresses. The vertical stresses Tyy used in the 
calculation include the stresses from 800 km seaward to 250 km landward of the 
trench. Assuming v — 0.25, we obtained the predicted topography shown in 
Table 2. Table 2 shows that for E = 70 GPa, which is a reasonable value for an 
oceanic lithosphere, the flexural parameter a is much smaller than the 
wavelength of the load applied at the base of the elastic lithosphere by viscous 
stresses, and the deflections due to the normal force are transmitted through 
the the elastic lithosphere essentially unaltered. The predicted topography of 
the elastic lithosphere due to the vertical stress Tyy does not agree with actual 
topography. 
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If the upper surface of the elastic lithosphere Is traction free, a stress 
due to the shear stress applied at the base of the elastic lithosphere by the 
viscous flow, must appear in the elastic lithosphere: 

= j£-/„ 0 rr( x ') dx ’ (36) 

where o%y is the shear stress on the base of the elastic lithosphere (Melosh and 
Raefsky 1980, equation (3)). From the viscous stresses obtained for model A, the 
calculated elastic stress (36) has a sign which is inconsistent with normal fault- 
ing in the outer trench region. 

Shown in Figure 12a, 12b, and 12c are contour plots for T**. Tyy and pres- 
sure, respectively. Figure 12d shows the orientation and magnitude of the prin- 
cipal compressional stresses in the slab at the center of each element. The 
stress pattern is very different from that of the circular arc slab model. From 
the top surface to the bottom surface of the lithosphere, compressional stress 
T*, increases in magnitude on the seaward side, while the pressure P and the 
vertical stress Tyy decrease. On the landward side, pressure is compressional on 
the top surface and tensional on the bottom surface. 

Another model with a different shape of the top surface but the same boun- 
dary conditions (model B) was tested. In model B, the top surface has a small 
dip angle. The difference between the two models is very small. Models with 
buoyancy force due to average density differences between the subducting slab 
and asthenosphere with a Newtonian viscosity were also considered. We assume 
there is a 0.0665 g/ cm 2 average density contrast between the cooling slab and 
the mantle (Hager Sc O’Connell 1981). The results near the trench are similar 
with those described above, and become more tensional at depth. 
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Observed grid with non-Newtonian Viscosity 

We also investigate the effects of the possible nonlinear dependence of 
viscosity on stresses. A constitutive law for non-Newtonian viscosity often used 
in the literature is: 


®ij I steady state = ( r, “ ) ff ij 


2 W 


(37) 


where is the strain rate, Tj e n is effective viscosity, which is stress dependent, 
and ay is the stress tensor. The effective viscosity for a power law fluid is 


V efl = f] n / a 


n-1 


(38) 


where a is the square root of the second deviatoric stress invariant and rj n is a 
constant (when n = 1, it is identical to the viscosity). Figure 13 shows the 
stresses for model A for a power law rheology with n = 3. The constant rj 3 does 
not depend upon depth. The general features of these plots are much different 
than those found for a Newtonian viscosity. In the region with small stresses 
seaward of the trench, the stress pattern becomes complicated, more asym- 
metric, and more sensitive to the velocity boundary conditions. The pressure 
near the trench becomes compressional and the region of tensional vertical 
stress Tyy moves farther away from the trench. This model does not predict the 
actual topography. 

By changing boundary conditions, we were able to find a model which 
predicts a vertical stress pattern similar to the trench topography. In this 
model, the slab begins bending on the seaward side of the trench, inconsistent 
with the actual topography. Since a very small change of the geometry of the 
model leads to significant variation of the stresses obtained, it is most unlikely 
that such a model is reasonable for the unique feature of trench topography of 


the world. 
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Con eluding Remarks 

The most important result of this work is that the analytic and numerical 
calculations for the stresses in a bending viscous flow show that the stress pat- 
tern is very sensitive to the velocity boundary conditions. The vertical stresses 
Tyy for non-Newtonian flow are more complex and asymmetric than for the 
Newtonian flow model. Regions where stresses are small are relatively more per- 
turbed than regions where stresses are large, and the transition from vertical 
tension to compression does not move in the direction to predict actual topogra- 
phy. 

For a non-Newtonian fluid with a power law 3 rheology, small stresses 
correspond to a higher viscosity, and are more sensitive to strain rate changes 
than for the Newtonian fluid. In our model, for the non-Newtonian flow, stresses 
near the trench are lower than for the deep part of the slab, and more complex 
than in Newtonian flow. This shows that in the region near the trench, the strain 
rate changes in the subducting slab, which are probably larger than in the other 
part of the slab, do not induce large stress variations as is expected in viscous 
flow. This suggests that topography near trenches is a solid deformation 
phenomenon rather than viscous flow Phenomenon. 

The primary purpose of this paper is to see if the first order topography and 
gravity field of subduction zones can be explained by viscous stresses generated 
in the lower viscous lithosphere during its subduction. Finite element models 
have been developed previously to explain the topographic structure of subduc- 
tion zones (Melosb & Raefsky 1980). However, the two models did not use realis- 
tic geometries. In their first model, the geometry of the upper surface of the 
viscous lithosphere is prescribed by error function profiles, which are not 
representative of actual subducting slabs. Although the stress Tyy obtained for 
the model is similar to observed topography, the sensitivity of stress distribu- 
tion to velocity boundary conditions, which are determined by the geometry of 
the slab, makes the result for any model with particular geometry not 
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applicable to other models, even approximately. Second, the results for the 
model with a circular arc velocity profile are complicated fay the smoothing pro- 
cedure, as explained in Part 1 of this paper. 

Summarizing the finite element computations for observed grids, we find: 
1) The stresses near the top surface are very sensitive to the velocity boundary 
conditions. 2) The viscous stresses do not have a pattern which is related even 
qualitatively to the topography near trenches, if the observed geometry is used 
to derive velocities. 

Similar problems caused by the sensitivity of stresses to the velocity boun- 
dary conditions exist in modeling the flow in the mantle outside the cold slab. 
Several authors used the comer flow model to investigate the torque balance on 
the slab (Stevenson and Turner 1977), mantle flow pressure and viscous coupling 
between the overriding plate and the down going slab (Tovish et al. 1978), and 
the compensation of geoid anomalies due to subducting slabs (McAdoo, 1982). 
When McKenzie (1969) first proposed to use corner flow to study the flow within 
the mantle, he emphasized that there was some doubt whether the resulting 
solutions apply to the flow in the earth, and in this respect there was a consider- 
able contrast between the section of the flow in the mantle and the previous sec- 
tions of his paper. In the corner flow model, there is a velocity jump at the 
trench which makes the solution singular at the trench axis. This singularity 
determines the whole solution pattern, and makes the solution inapplicable 
quantitatively to actual subduction zones in the distance range of the length of 
the slab, although it is valid qualitatively. 

The results of our models using observed geometries as constraints show 
diverse stress patterns, out of phase with the actual topography and gravity. 
The conclusion is that the topography and geoid anomalies associated with 
trenches are supported by the strength of the strong upper lithosphere which 
behaves as an elastic-plastic solid at the appropriate time scale. 
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Topography appears to be an effective discriminant in determining the 
appropriate rheological model for the lithosphere. Study of the topography of 
an actual subduction zone suggests an elastic-plastic rheology for oceanic litho- 
sphere. More precise results may be obtained by applying the finite element 
method to an elastic-plastic layer overlaying a viscous layer. It appears that 
sophisticated methods and accurate data are necessary for a better under- 
standing of the lithosphere’s behavior. 

Although the sensitivity of stress to the velocity boundary condition dimin- 
ishes the possibility that viscous stress supports the topography, the viscous 
behavior of subducting slab is obvious in many aspects as discussed in the intro- 
duction. Viscous stresses acting on the elastic lithosphere may play a more 
important role in the dipping section than in the bending and horizontal section 
of the elastic lithosphere. 

Consider the lithosphere model described in Figure la of Melosh 1978. In 
our calculations, we find that the stress' pattern in the viscous slab after bending 
is stable. In Figure 12, the deviatoric stresses show that the tensional axis is 
parallel to the slab. This result is in accord with the existence of tensional 
events in the lower part of the upper lithosphere. Small changes of velocities 
caused by the interaction of two converging plates may enhance or diminish the 
tensional stress. Enhanced tensional stresses in the lower surface of the upper 
lithosphere will make upper compressional events become tensional or make 
the double seismic zone a tensional zone. Diminished tensional stresses in the 
lower surface of the upper lithosphere are favorable a single one compressional 
zone. 
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Figure captions 


Figure 1. (a) Geometry for corner flow in region A. (b) The geometry of a 
simple subducting slab model with dip angle 6 . The slab has a constant 
thickness in the direction normal to the top surface, which is composed 
of two straight lines connecting a circular arc DC with radius b and center 
A, The bottom surface is concentric with the top surface and has radius 
a. 


Figure 2. Pressure profiles measured at two different dimensionless depths d 
in the bending slab model. Dashed lines show p f , the pressure which is 
calculated by finite element methods and which has a minimum value ( 
tensional ) at x = 0 ( Fig 1(b) ). Solid lines show p a , the pressure which is 
calculated with the analytic formula (32). The scales are different for 
each depth, and are normalized relative to the minimum pressure. 

Figure 3. (a) Pressure contours for the circular-arc grid with constant 

tangential velocity boundary condition on the top surface. 1470 iso- 
parametric bilinear elements, distributed according to the numbers on 
the sides of the grid, were used. The velocity vectors on the beginning 
and end of the grid are parallel to the upper surface of the slab and have 
unit magnitude. The bottom surface of the slab is traction free, (b) Pres- 
sure profile at the upper surface of the slab. 

Figure 4. (a) Pressure contours for a circular arc grid with smoothed velocity 
boundary conditions on top surface. The bottom surface boundary condi- 
tions are the same as those stated in Figure 4. (b) Pressure profile for 
the layer at distance d = 0.005 from the upper surface. 


Figure 5. Three error function velocity profiles (a-1), (a-2), and (a-3), 
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determined from formula (35) with v = 1.0, 6 = 45°, and 

S c = 0.8, 0.4, and 0.2 respectively, (b-l), (b-2), and (b-3) are the pressure 
profiles along the top line of the circular-arc grid with velocity profiles 
(a-1), (a-2), and (a-3) applied respectively, (c-l), (c-2), and (c-3) are for 
the rectangular grid with the same velocity boundary conditions (a-l), 
(a-2), and (a-3). 

Figure 6. The velocity vectors of a small segment of the stream line pro- 
jected onto two characteristic directions. Near the bend axis the tangen- 
tial velocity gradient is zero and the normal velocity gradient is constant. 

Figure 7. Superposition of the pressure due to normal gradients. The 10 thin 
lines are pressure profiles calculated using the fundamental solution in 
Table 1 for a vertical step function velocity boundary condition. The 
parameters used in the formula are 77 = 1.0, y = -0.005, and 
X = 0, ±1.0, ±0.75, ±0.50, ±0.25. The thick line is a pressure profile which 
sums these 10 individual pressure profiles. 

Figure 8. The pressure along the top line of the circular arc slab model cal- 
culated using the finite element method, (a - 1) Velocity boundary condi- 
tions as a function of horizontal distance applied to the circular arc boun- 
dary of the grid shown as Figure 3(a), where the thickness H v = 1, radius 
R = 1.85, dip angle 6 = 45°, and where ( x = 0 ) corresponds to the begin- 
ning of the arc. A, and Ay are the horizontal ( seaward ) and vertical ( 
upward ) velocity components, v x and v y in formula (33) of the text, 
respectively. These velocity profiles have discontinuous first derivatives 
at the two ends of the circular arc section. Curves B x and By are 
smoothed versions of A x and A y ( i.e. v x and v y in formula (34) of the 
text, which have a continuous first derivative ). C x and C y are the 
differences between B x and A x , and B y and A y respectively, (a-2) The 
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differential velocities C x - and C y - plotted on a coordinate system which has 
the origin at the middle of the circular arc. (b - 1), (b - 2), (b - 3) are 
pressure profiles along the top surface of the circular arc grid, which 
correspond to the velocity boundary conditions A* and Ay, C x and C y , and 
B x and B y respectively, (c - 1) and (d - 1) are schematic normal velocity 
gradients gj. for the velocity boundary condition of A* and Ay. and g' x for 
the differential velocity boundary condition C x and C y . The pressure 
profiles for the normal velocity gradients and g' a are shown in (c - 2) 
and (d - 2). (e) A diagram which illustrates how to obtain the normal velo- 
city gradient gj_ and tangential velocity gradient gj | for a constant tangen- 
tial velocity boundary condition on the circular arc. The tangential velo- 
city gradient g|| is zero, and the normal velocity gradient g^is a constant 
on the the entire arc. 

Figure 9. Two finite element grids used in the present calculations, which 
have slightly different geometries for the top line of the grid. The hor- 
izontal solid lines above each grid are the x-axis (at sea level). The trench 
axis off Northern Honshu, Japan is located at x = 0. The dotted line 
beneath the x-axis is the averaged topography data of Northern Honshu. 
Hypocenters of low-angle thrust and down-dip compressional events of 
Northern Honshu, Japan are plotted as dots below the topography line. 
Above the grids the horizontal velocities v x (landward) and vertical veloci- 
ties v y (downward) are plotted. The top line of the grid in (a), called 
model A in the text, is a polynomial obtained by a least squares fit to the 
topography and upper envelope of the hypocenters. The top line of the 
grid in (b), called model B. is a polynomial which fits the topography and 
upper envelope of the hypocenters slightly differently. 

Figure 10. The bathymetric profile east of Northern Honshu. The seaward 
side of the profile is obtained by averaging 5 profiles ranging from 450 km 
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west to 900 km east of the trench axis. These cross sections correspond 
to the boundaries of the 4 seismicity cross sections used by Kawakatsu 
and Seno (1983). The top lines of the two finite element grids, for model A 
and model B are also plotted. 

Figure 11. The stress components T Ml Tyy, and pressure P as functions of dis- 
tance from trench for the top (a) and bottom (b) surface of the model A 
with a Newtonian constitutive law, 

Figure 12. The contours for stress components T D (a), Tyy (b), and pressure P 
(c). and the deviatoric compressional principal stress (d) for model A with 
Newtonian constitutive law. The length of the small lines plotted in (d) is 
proportional to the magnitude of the compressional principal stress and 
their orientation is parallel to the direction of the compressional princi- 
pal stress. 

Figure 13. The stress components T^.Tyy, and pressure P as functions of dis- 
tance from trench for the top (a) and bottom (b) surface of the model A 
with a non-Newtonian constitutive law ( power law 3 ). 



Table 1 

Fundamental Solutions 
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Boundary condition 
v t (x,0) = CdsXjc, 
v y (x,0) = 0, A>0 

Stream function +(x,y) 
ye Xy Cos Ax 

Pressure p{x ,y)-p{x ,-») 
2T)Xe * v SinAx 


J>(*.v)-p (*.—■) 
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Boundary condition 
v, (x .0) = 0. 

Vy(x.O) = Cbs\x. A>0 

Stream function +(x.y) 
(y-^-)e Xv SinAx 

Pressure p{x .y)-p{x .-») 
—2r)\e Xv Cos Ax 


P(*.V)“P (*.-“) 
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Boundary condition 
v,(x.O) = «(x-x 0 ). 
v„(x.0) = 0 

Stream function +(x.y) 


n (x-x 0 ) 2 +y 2 

Pressure j>(x ,y)-p{x ,-») 
_jl 4(z-*o)v 
w [(*~*o) 2+ y*] 2 


p(*.v)-p(*.— ) 
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Boundary condition 
v t (x,0)=0. 
v„(x.O)=<5(x-x 0 ) 


Stream function +(x,y) 
1 f (I-I 0 )v 


- I(x-x 0 ) 2 +y 2 


— arccot 


*~*o 

V 


Pressu re p (x .y)-p (x ,— ) 
2r? y g -(x -x 0 ) 2 
" tt [(x-x 0 ) 2 + y 2 } 2 
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Boundary condition 
v I (x,0) = H (x-x 0 ). 
v y (x ,0)=0 

Stream function +(x.y) 

1 ,*-*o 

— y arccot 

n y 

Pressure p(x .y)-p(x,-=>) 
ZL. 2v - 
n (x-x 0 ) 2 + V 2 
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Boundary condition 
v,(x.0) = 0. 

Vy(x.O) = H(x-x 0 ) 

Stream function +(x,y) 
*-*o 




■x 0 ) arccof- 


Pressure p (x ,y)-p (x 

T 7 -»p) 

~n (x-x 0 ) 2 + y 2 


P .-**) 


p(*.v)-p(* .-•) 



p(*.v)-p (*.-«■) 





Table 2 



x t km 
w t m 
x r km 
w r m 


Predicted Topography 



Elastic plate 



over viscous plate 
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Viscous 
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E = 0.7x10” Pa 


plate 

H e = 15 km 
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topography 
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D c =1.7xl0 23 
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